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Abstract 

We present evidenee that the prismatie and seeondary prism faeets of iee-Ih erystals possess struetural 
features that ean reduee the effeetive hydrophilieity of the iee/water interfaee. The spreading dynamies of 
liquid water droplets on iee faeets exhibits long-time behavior that differs for the prismatie {1010} and 
seeondary prism {1120} faeets when eompared with the basal {0001} and pyramidal {2021} faeets. We 
also present the results of simulations of solid-liquid frietion of the same four erystal faeets being drawn 
through liquid water, and find that the two prismatie faeets exhibit roughly half the solid-liquid frietion 
of the basal and pyramidal faeets. These simulations provide evidenee that the two prismatie faees have a 
signifieantly smaller effeetive surfaee area in eontaet with the liquid water. The iee / water interfaeial widths 
for all four erystal faeets are similar (using both struetural and dynamie measures), and were found to be 
independent of the shear rate. Additionally, deeomposition of orientational time eorrelation funetions show 
position-dependenee for the short- and longer-time deeay eomponents elose to the interfaee. 
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I. INTRODUCTION 


Surfaces can be characterized as hydrophobic or hydrophilic based on the strength of the in¬ 
teractions with water. Hydrophobic surfaces do not have strong enough interactions with water 
to overcome the internal attraction between molecules in the liquid phase, and the degree of hy- 
drophilicity of a surface can be described by the extent a droplet can spread out over the surface. 
The contact angle, 0, formed between the solid and the liquid depends on the free energies of the 
three interfaces involved, and is given by Young’s equation [1], 

cos0 = (y,v-Y^/)/Y/v (1) 

Here Jsv, Jsi, and Yzv are the free energies of the solid/vapor, solid/liquid, and liquid/vapor inter¬ 
faces, respectively. Large contact angles, 0 > 90°, correspond to hydrophobic surfaces with low 
wettability, while small contact angles, 0 < 90°, correspond to hydrophilic surfaces. Experimen¬ 
tally, measurements of the contact angle of sessile drops is often used to quantify the extent of 
wetting on surfaces with thermally selective wetting characteristics [2-4]. 

Nanometer-scale structural features of a solid surface can influence the hydrophilicity to a 
surprising degree. Small changes in the heights and widths of nano-pillars can change a surface 
from superhydrophobic, 0 > 150°, to hydrophilic, 0 ~ 0° [5]. This is often referred to as the 
Cassie-Baxter to Wenzel transition. Nano-pillared surfaces with electrically tunable Cassie-Baxter 
and Wenzel states have also been observed [6-10]. Luzar and coworkers have modeled these 
transitions on nano-patterned surfaces [11-14], and have found the change in contact angle is due 
to the field-induced perturbation of hydrogen bonding at the liquid/vapor interface [11]. 

One would expect the interfaces of ice to be highly hydrophilic (and possibly the most hy¬ 
drophilic of all solid surfaces). In this paper we present evidence that some of the crystal facets 
of ice-Ih have structural features that can reduce the effective hydrophilicity. Our evidence for this 
comes from molecular dynamics (MD) simulations of the spreading dynamics of liquid droplets 
on these facets, as well as reverse non-equilibrium molecular dynamics (RNEMD) simulations of 
solid-liquid friction. 

Quiescent ice-Ih/water interfaces have been studied extensively using computer simulations. 
Hayward and Haymet characterized and measured the widths of these interfaces [15, 16]. Nada 
and Eurukawa have also modeled the width of basal/water and prismatic/water interfaces [17] as 
well as crystal restructuring at temperatures approaching the melting point [18]. 
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The surface of ice exhibits a pre-melting layer, often called a quasi-liquid layer (QLL), at 
temperatures near the melting point. MD simulations of the facets of ice-Ih exposed to vacuum 
have found QLL widths of approximately 10 A at 3 K below the melting point [19]. Similarly, 
Limmer and Chandler have used the mW water model [20] and statistical field theory to estimate 
QLL widths at similar temperatures to be about 3 nm [21]. 

Recently, Sazaki and Furukawa have developed a technique using laser confocal microscopy 
combined with differential interference contrast microscopy that has sufficient spatial and tempo¬ 
ral resolution to visualize and quantitatively analyze QLLs on ice crystals at temperatures near 
melting [22]. They have found the width of the QLLs perpendicular to the surface at -2.2"C to be 
3-4 A wide. They have also seen the formation of two immiscible QLLs, which displayed different 
dynamics on the crystal surface [23]. 

Using molecular dynamics simulations, Samadashvili has recently shown that when two 
smooth ice slabs slide past one another, a stable liquid-like layer develops between them [24]. 
In a previous study, our RNEMD simulations of ice-Ih shearing through liquid water have pro¬ 
vided quantitative estimates of the solid-liquid kinetic friction coefficients [25]. These displayed a 
factor of two difference between the basal and prismatic facets. The friction was found to be inde¬ 
pendent of shear direction relative to the surface orientation. We attributed facet-based difference 
in liquid-solid friction to the 6.5 A corrugation of the prismatic face which reduces the effective 
surface area of the ice that is in direct contact with liquid water. 

In the sections that follow, we describe the simulations of droplet-spreading dynamics using 
standard MD as well as simulations of tribological properties using RNEMD. These simulations 
give complementary results that point to the prismatic and secondary prism facets having roughly 
half of their surface area in direct contact with the liquid. 


IL METHODOLOGY 

A, Construction of the Ice / Water interfaces 

Ice Ih crystallizes in the hexagonal space group P 63 /mmc, and common ice crystals form hexag¬ 
onal plates with the basal face, { 0001 }, forming the top and bottom of each plate, and the prismatic 
facet, {1010}, forming the sides. In extreme temperatures or low water saturation conditions, ice 
crystals can easily form as hollow columns, needles and dendrites. These are structures that ex- 
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pose other erystalline faeets of the iee to the surroundings. Among the more eommon faeets are 
the seeondary prism, {1120}, and pyramidal, {2021}, faees. 

We found it most useful to work with proton-ordered, zero-dipole erystals that expose strips 
of dangling H-atoms and lone pairs [26]. Our struetures were ereated starting from Strueture 6 of 
Hirseh and Ojamae’s set of orthorhombie representations for iee-I/, [27]. This erystal strueture was 
eleaved along the four different faeets. The exposed faee was reoriented normal to the z-axis of the 
simulation eell, and the struetures were extended to form large exposed faeets in reetangular box 
geometries. Liquid water boxes were ereated with identieal dimensions (in x and y) as the iee, with 
a z dimension of three times that of the iee bloek, and a density eorresponding to 1 g / em^. Eaeh 
of the iee slabs and water boxes were independently equilibrated at a pressure of 1 atm, and the 
resulting systems were merged by earving out any liquid water moleeules within 3 A of any atoms 
in the iee slabs. Eaeh of the eombined iee/water systems were then equilibrated at 225K, whieh 
is the liquid-iee eoexistenee temperature for SPC/E water [28]. Referenee 25 eontains a more 
detailed explanation of the eonstruetion of similar iee/water interfaees. The resulting dimensions 
as well as the number of iee and liquid water moleeules eontained in eaeh of these systems are 
shown in Table I. 

The SPC/E water model [29] has been extensively eharaeterized over a wide range of liquid 
eonditions [30, 31], and its phase diagram has been well studied [32-35]. With longer eutoff radii 
and eareful treatment of eleetrostaties, SPC/E mostly avoids metastable erystalline morphologies 
like iee-/ [35] and iee-B [32]. The free energies and melting points [28, 30, 32-39] of various 
other erystalline polymorphs have also been ealeulated. Haymet et al. have studied quieseent 
lee-Ih/water interfaees using the SPC/E water model, and have seen struetural and dynamie mea¬ 
surements of the interfaeial width that agree well with more expensive water models, although 
the eoexistenee temperature for SPC/E is still well below the experimental melting point of real 
water [28]. Given the extensive data and speed of this model, it is a reasonable ehoiee even though 
the temperatures required are somewhat lower than real iee / water interfaees. 


Ill, DROPLET SIMULATIONS 

lee surfaees with a thiekness of ~ 20 A were ereated as deseribed above, but were not sol¬ 
vated in a liquid box. The erystals were then replieated along the x and y axes (parallel to the 
surfaee) until a large surfaee (> 126 nm^) had been ereated. The sizes and numbers of moleeules 
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in each of the surfaces is given in Table SI. Weak translational restraining potentials with spring 
constants of 1.5 kcal moP A (prismatic and pyramidal facets) or 4.0 kcal moP A (basal 
facet) were applied to the centers of mass of each molecule in order to prevent surface melting, 
although the molecules were allowed to reorient freely. A water droplet containing 2048 SPC/E 
molecules was created separately. Droplets of this size can produce agreement with the Young 
contact angle extrapolated to an infinite drop size [12]. The surfaces and droplet were indepen¬ 
dently equilibrated to 225 K, at which time the droplet was placed 3-5 A above the surface. Five 
statistically independent simulations were carried out for each facet, and the droplet was placed at 
unique x and y locations for each of these simulations. Each simulation was 5 ns in length and was 
conducted in the microcanonical (NVE) ensemble. Representative configurations for the droplet 
on the prismatic facet are shown in figure 1. 

IV. SHEARING SIMULATIONS (INTERFACES IN BULK WATER) 

To perform the shearing simulations, the velocity shearing and scaling variant of reverse non¬ 
equilibrium molecular dynamics (VSS-RNEMD) was employed [31]. This method performs a se¬ 
ries of simultaneous non-equilibrium exchanges of linear momentum and kinetic energy between 
two physically-separated regions of the simulation cell. The system responds to this unphysical 
flux with velocity and temperature gradients. When VSS-RNEMD is applied to bulk liquids, trans¬ 
port properties like the thermal conductivity and the shear viscosity are easily extracted assuming 
a linear response between the flux and the gradient. At the interfaces between dissimilar materials, 
the same method can be used to extract interfacial transport properties (e.g. the interfacial thermal 
conductance and the hydrodynamic slip length). 

The kinetic energy flux (producing a thermal gradient) is necessary when performing shearing 
simulations at the ice-water interface in order to prevent the frictional heating due to the shear from 
melting the crystal. Reference 25 provides more details on the VSS-RNEMD method as applied to 
ice-water interfaces. A representative configuration of the solvated prismatic facet being sheared 
through liquid water is shown in figure 2. 

All simulations were performed using OpenMD [40, 41], with a time step of 2 fs and periodic 
boundary conditions in all three dimensions. Electrostatics were handled using the damped-shifted 
force real-space electrostatic kernel [42]. 

The interfaces were equilibrated for a total of 10 ns at equilibrium conditions before being 
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exposed to either a shear or thermal gradient. This eonsisted of 5 ns under a eonstant temperature 
(NVT) integrator set to 225 K followed by 5 ns under a mieroeanonieal (NVE) integrator. Weak 
thermal gradients were allowed to develop using the VSS-RNEMD (NVE) integrator using a small 
thermal flux (—2.0 x 10^^ keal/mol/A^/fs) for a duration of 5 ns to allow the gradient to stabilize. 
The resulting temperature gradient was ~ lOK over the entire box length, whieh was suffieient to 
keep the temperature at the interfaee within ± 1 K of the 225 K target. 

Veloeity gradients were then imposed using the VSS-RNEMD (NVE) integrator with a range of 
momentum fluxes. The systems were divided into 100 bins along the z-axis for the VSS-RNEMD 
moves, whieh were attempted every time step. Although eomputationally expensive, this was 
done to minimize the magnitude of eaeh individual momentum exehange. Beeause individual 
VSS-RNEMD exehange moves eonserve both total energy and linear momentum, the method ean 
be “bolted-on” to simulations in any ensemble. The simulations of the pyramidal interfaee were 
performed under the eanonieal (NVT) ensemble. When time eorrelation funetions were eomputed, 
the RNEMD simulations were done in the mieroeanonieal (NVE) ensemble. All simulations of 
the other interfaees were earried out in the mieroeanonieal ensemble. 

These gradients were allowed to stabilize for 1 ns before data eolleetion started. Onee estab¬ 
lished, four sueeessive 0.5 ns runs were performed for eaeh shear rate. During these simulations, 
eonfigurations of the system were stored every 1 ps, and statisties on the strueture and dynamies 
in eaeh bin were aeeumulated throughout the simulations. Although there was some small varia¬ 
tion in the measured interfaeial width between sueeeessive runs, no indieation of bulk melting or 
erystallization was observed. That is, no large seale ehanges in the positions of the top and bottom 
interfaees oeeurred during the simulations. 


V. RESULTS 

A, Ice - Water Contact Angles 

To determine the extent of wetting for eaeh of the four erystal faeets, eontaet angles for liquid 
droplets on the iee surfaees were eomputed using two methods. In the first method, the droplet is 
assumed to form a spherieal eap, and the eontaet angle is estimated from the z-axis loeation of the 
droplet’s eenter of mass (zcm)- This proeedure was first deseribed by Hautman and Klein [43], and 
was utilized by Hirvi and Pakkanen in their investigation of water droplets on polyethylene and 


6 



poly(vinyl chloride) surfaces [44]. For each stored configuration, the contact angle, 0, was found 
by inverting the expression for the location of the droplet center of mass. 


(Zcm) = 


1 — cos9 3 + cos9 

2 + cos9 J 2 + cos9 ’ 


( 2 ) 


where Rq is the radius of the free water droplet. 

In addition to the spherical cap method outlined above, a second method for obtaining the 
contact angle was described by Ruijter, Blake, and Coninck [45]. This method uses a cylindrical 
averaging of the droplet’s density profile. A threshold density of 0.5 g cm'^ is used to estimate 
the location of the edge of the droplet. The r and z-dependence of the droplet’s edge is then fit 
to a circle, and the contact angle is computed from the intersection of the fit circle with the z-axis 
location of the solid surface. Again, for each stored configuration, the density profile in a set of 
annular shells was computed. Due to large density fluctuations close to the ice, all shells located 
within 2 A of the ice surface were left out of the circular fits. The height of the solid surface 
(zsuface) along with the best fitting origin (zdropiet) and radius (rdropiet) of the droplet can then be 
used to compute the contact angle, 

f. nn I • -1 ( 2:droplet ~^surface \ 

0 = 90 H-sm - - - . (3) 

tt \ ^droplet / 

Both methods provided similar estimates of the dynamic contact angle, although the spherical cap 
method is significantly less prone to noise, and is the method used to compute the contact angles 
in table II. 

Because the initial droplet was placed above the surface, the initial value of 180° decayed over 
time (See figure 1 in the SI). Each of these profiles were fit to a biexponential decay, with a short- 
time contribution (x^) that describes the initial contact with the surface, a long time contribution 
(Xj) that describes the spread of the droplet over the surface, and a constant (Ooo) to capture the 
infinite-time estimate of the equilibrium contact angle. 


0(t) = 000 +(180-Ooo) 


ae A'^" + (l-a)e 


(4) 


We have found that the rate for water droplet spreading across all four crystal facets, ^spread = 
l/x^ ~ 0.7 ns^^ However, the basal and pyramidal facets produced estimated equilibrium contact 
angles, 0oo ~ 35°, while prismatic and secondary prismatic had values for 0oo near 43° as seen in 
Table II. 
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These results indicate that by traditional measures, the basal and pyramidal facets are more 
hydrophilic than the prismatic and secondary prism facets, and surprisingly, that the differential 
hydrophilicities of the crystal facets is not reflected in the spreading rate of the droplet. 

B, Solid-liquid friction of the interfaces 

In a bulk fluid, the shear viscosity, T], can be determined assuming a linear response to a shear 
stress, 

jz{Px)=V^. (5) 

Here jz{Px) is the flux (in x-momentum) that is transferred in the z direction (i.e. the shear stress). 
The RNEMD simulations impose an artificial momentum flux between two regions of the simula¬ 
tion, and the velocity gradient is the fluid’s response. This technique has now been applied quite 
widely to determine the viscosities of a number of bulk fluids [46-48]. 

At the interface between two phases (e.g. liquid / solid) the same momentum flux creates a 
velocity difference between the two materials, and this can be used to define an interfacial friction 
coefficient (k), 

jziPx) = K [Vxiliquid) - Vxisolid)] (6) 

where Vx{solid) is the velocity of the solid and Vx{liquid) is the velocity of the liquid measured at 
the hydrodynamic boundary layer. 

The simulations described here contain significant quantities of both liquid and solid phases, 
and the momentum flux must traverse a region of the liquid that is simultaneously under a thermal 
gradient. Since the liquid has a temperature-dependent shear viscosity, ri(r), estimates of the 
solid-liquid friction coefficient can be obtained if one knows the viscosity of the liquid at the 
interface (i.e. at the melting temperature, Tm), 

^ ^ _ f .rj. 

[vxifluid) - Vx{solid)] \^z J ' 

For SPC/E, the melting temperature of Ice-Ih is estimated to be 225 K [28]. To obtain the value 
of ri(225 K) for the SPC/E model, a 31.09 x 29.38 x 124.39 A box with 3744 water molecules in 
a disordered configuration was equilibrated to 225 K, and five statistically-independent shearing 
simulations were performed (with imposed fluxes that spanned a range of 3 —)■ 13 m s~^ ). Each 
simulation was conducted in the microcanonical ensemble with total simulation times of 5 ns. 
The VSS-RNEMD exchanges were carried out every 2 fs. We estimate ri(225 K) to be 0.0148 ± 
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0.0007 Pa s for SPC/E, roughly ten times larger than the shear viscosity previously computed at 
280 K [31], 

The interfacial friction coefficient can equivalently be expressed as the ratio of the viscosity of 
the fluid to the hydrodynamic slip length, K = ri/5. The slip length is an indication of strength of 
the interactions between the solid and liquid phases, although the connection between slip length 
and surface hydrophobicity is not yet clear. In some simulations, the slip length has been found 
to have a link to the effective surface hydrophobicity [49], although Ho et al. have found that 
liquid water can also slip on hydrophilic surfaces [50]. Experimental evidence for a direct tie be¬ 
tween slip length and hydrophobicity is also not definitive. Total-internal reflection particle image 
velocimetry (TIR-PIV) studies have suggested that there is a link between slip length and effec¬ 
tive hydrophobicity [51, 52]. However, recent surface sensitive cross-correlation spectroscopy 
(TIR-ECCS) measurements have seen similar slip behavior for both hydrophobic and hydrophilic 
surfaces [53]. 

In each of the systems studied here, the interfacial temperature was kept fixed to 225 K, which 
ensured the viscosity of the fluid at the interace was identical. Thus, any significant variation in 
K between the systems is a direct indicator of the slip length and the effective interaction strength 
between the solid and liquid phases. 

The calculated k values found for the four crystal facets of Ice-Ih are shown in Table II. The 
basal and pyramidal facets were found to have similar values of k ^ 6 (x lO^^amu A ^fs~^), while 
the prismatic and secondary prism facets exhibited K ^ 3 (x lO^^amu A ^fs^^). These results are 
also essentially independent of the direction of the shear relative to channels on the surfaces of 
the facets. The friction coefficients indicate that the basal and pyramidal facets have significantly 
stronger interactions with liquid water than either of the two prismatic facets. This is in agreement 
with the contact angle results above - both of the high-friction facets exhibited smaller contact 
angles, suggesting that the solid-liquid friction (and inverse slip length) is correlated with the 
hydrophilicity of these facets. 

C. Structural measures of interfacial width under shear 

One of the open questions about ice/water interfaces is whether the thickness of the ’slush¬ 
like’ quasi-liquid layer (QEE) depends on the facet of ice presented to the water. In the QEE 
region, the water molecules are ordered differently than in either the solid or liquid phases, and 
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also exhibit distinct dynamical behavior. The width of this quasi-liquid layer has been estimated 
by finding the distance over which structural order parameters or dynamic properties change from 
their bulk liquid values to those of the solid ice. The properties used to find interfacial widths have 
included the local density, the diffusion constant, and both translational and orientational order 
parameters [15, 16, 25, 28, 36, 54, 55]. 

The VSS-RNEMD simulations impose thermal and velocity gradients. These gradients per¬ 
turb the momenta of the water molecules, so parameters that depend on translational motion are 
often measuring the momentum exchange, and not physical properties of the interface. As a struc¬ 
tural measure of the interface, we have used the local tetrahedral order parameter, which mea¬ 
sures the match of the local molecular environments (e.g. the angles between nearest neighbor 
molecules) to perfect tetrahedral ordering. This quantity was originally described by Errington 
and Debenedetti [56] and has been used in bulk simulations by Kumar et al. [57] It has previously 
been used in ice/water interfaces by by Bryk and Haymet [33]. 

To determine the structural widths of the interfaces under shear, each of the systems was di¬ 
vided into 100 bins along the z-dimension, and the local tetrahedral order parameter (Eq. 5 in 
Reference 25) was time-averaged in each bin for the duration of the shearing simulation. The 
spatial dependence of this order parameter, q{z), is the tetrahedrality profile of the interface. The 
lower panels in figures 2-5 in the supporting information show tetrahedrality profiles (in circles) 
for each of the four interfaces. The q{z) function has a range of (0,1), where a value of unity 
indicates a perfectly tetrahedral environment. The q{z) for the bulk liquid was found to be ~ 0.77, 
while values of ~ 0.92 were more common in the ice. The tetrahedrality profiles were fit using 
a hyperbolic tangent function (see Eq. 6 in Reference 25) designed to smoothly fit the bulk to ice 
transition while accounting for the weak thermal gradient. In panels b and c of the same figures, 
the resulting thermal and velocity gradients from an imposed kinetic energy and momentum fluxes 
can be seen. The vertical dotted lines traversing these figures indicate the midpoints of the inter¬ 
faces as determined by the tetrahedrality profiles. The hyperbolic tangent fit provides an estimate 
of ^struct, the structural width of the interface. 

We find the interfacial width to be 3.2 ± 0.2 A (pyramidal) and 3.2 ± 0.2 A (secondary prism) 
for the control systems with no applied momentum flux. This is similar to our previous results for 
the interfacial widths of the quiescent basal (3.2 ± 0.4 A) and prismatic systems (3.6 ± 0.2 A). 

Over the range of shear rates investigated, 0.4 —)■ 6.0 ms^^ for the pyramidal system and 
0.6 —)■ 5.2 m s^^ for the secondary prism, we found no significant change in the interfacial width. 
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The mean interfacial widths are collected in table 11. This follows our previous findings of the 
basal and prismatic systems, in which the interfacial widths of the basal and prismatic facets were 
also found to be insensitive to the shear rate [25]. 

The similarity of these interfacial width estimates indicate that the particular facet of the ex¬ 
posed ice crystal has little to no effect on how far into the bulk the ice-like structural ordering 
persists. Also, it appears that for the shearing rates imposed in this study, the interfacial width is 
not structurally modified by the movement of water over the ice. 

D. Dynamic measures of interfacial width under shear 

The spatially-resolved orientational time correlation function, 

Czfct) = {Piium ■uKt))5(zK0) -z)), (8) 

provides local information about the decorrelation of molecular orientations in time. Here, P 2 is the 
second-order Legendre polynomial, and u/ is the molecular vector that bisects the HOH angle of 
molecule i. The angle brackets indicate an average over all the water molecules and time origins, 
and the delta function restricts the average to specific regions. In the crystal, decay of C 2 (z,t) 
is incomplete, while in the liquid, correlation times are typically measured in ps. Observing the 
spatial-transition between the decay regimes can define a dynamic measure of the interfacial width. 

To determine the dynamic widths of the interfaces under shear, each of the systems was di¬ 
vided into bins along the z-dimension 3 A wide) and C 2 (z,t) was computed using only those 
molecules that were in the bin at the initial time. To compute these correlation functions, each 
of the 0.5 ns simulations was followed by a shorter 200 ps microcanonical (NVE) simulation in 
which the positions and orientations of every molecule in the system were recorded every 0.1 ps. 

The time-dependence was fit to a triexponential decay, with three time constants: Tshorts mea¬ 
suring the librational motion of the water molecules, ^middle, measuring the timescale for breaking 
and making of hydrogen bonds, and Xiong, corresponding to the translational motion of the water 
molecules. An additional constant was introduced in the fits to describe molecules in the crystal 
which do not experience long-time orientational decay. 

In Figures 6-9 in the supporting information, the z-coordinate profiles for the three decay con¬ 
stants, Xshort, Xmiddie, and Xiong for the different interfaces are shown. (Figures 6 & 7 are new 
results, and Figures 8 & 9 are updated plots from Ref 25.) In the liquid regions of all four inter- 
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faces, we observe Xmiddie and X/ong to have approximately eonsistent values of 3 — 6 ps and 30 — 40 
ps, respeetively. Both of these times inerease in value approaehing the interfaee. Approaehing 
the interfaee, we also observe that Xshort deereases from its liquid-state value of 72 — 76 fs. The 
approximate values for the deeay eonstants and the trends approaehing the interfaee mateh those 
reported previously for the basal and prismatie interfaees. 

We have estimated the dynamie interfaeial width by fitting the profiles of all the three 
orientational time eonstants with an exponential deeay to the bulk-liquid behavior, 

x{z) ^ Xliqiii(l(x^,all ^ wall)/ dyn 

where Xuquid and x^aii are the liquid and projeeted wall values of the deeay eonstants, Zwaii is the 
loeation of the interfaee, as measured by the struetural order parameter. These values are shown 
in table IT Beeause the bins must be quite wide to obtain reasonable profiles of C 2 (z, t), the error 
estimates for the dynamie widths of the interfaee are signifieantly larger than for the struetural 
widths. However, all four interfaees exhibit dynamie widths that are signifieantly below 1 nm, and 
are in reasonable agreement with the struetural width above. 


VI, CONCLUSIONS 

In this work, we used MD simulations to measure the advaneing eontaet angles of water 
droplets on the basal, prismatie, pyramidal, and seeondary prism faeets of lee-Ih. Although we 
saw no signifieant ehange in the rate at whieh the droplets spread over the surfaee, the long-time 
behavior prediets larger equilibrium eontaet angles for the two prismatie faeets. 

We have also used RNEMD simulations of water interfaees with the same four erystal faeets 
to eompute solid-liquid Motion ooeffioients. We have observed ooeffioients of frietion that differ 
by a faetor of two between the two prismatie faeets and the basal and pyramidal faeets. Beeause 
the solid-liquid Motion ooeffioient is direotly tied to the inverse of the hydrodynamio slip length, 
this suggests that there are signifieant differenees in the overall interaetion strengths between these 
faeets and the liquid layers immediately in eontaet with them. 

The agreement between these two measures have lead us to eonelude that the two prismatie 
faeets have a lower hydrophilieity than either the basal or pyramidal faeets. One possible expla¬ 
nation of this behavior is that the faoe presented by both prismatie faeets oonsists of deep, narrow 
ohannels (i.e. stripes of adjaeent rows of pairs of hydrogen-bound water moleoules). At the sur- 
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faces of these faeets, the ehannels are 6.35 A wide and the sub-surfaee iee layer is 2.25 A below 
(and therefore bloeked from hydrogen bonding with the liquid). This means that only 1/2 of the 
surfaee moleeules ean form hydrogen bonds with liquid-phase moleeules. 

In the basal plane, the surfaee features are shallower (1.3 A), with no bloeked subsurfaee layer. 
The pyramidal faee has mueh wider ehannels (8.65 A) whieh are also quite shallow (1.37 A). These 
features allow liquid phase moleeules to form hydrogen bonds with all of the surfaee moleeules 
in the basal and pyramidal faeets. This means that for similar surfaee areas, the two prismatie 
faeets have an effeetive hydrogen bonding surfaee area of half of the basal and pyramidal faeets. 
The reduetion in the effeetive surfaee area would explain mueh of the behavior observed in our 
simulations. 
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FIG. 1: Computational model of a droplet of liquid water spreading over the prismatic {lOlO} facet of ice, 
before (left) and 2.6 ns after (right) being introduced to the surface. The contact angle (0) shrinks as the 
simulation proceeds, and the long-time behavior of this angle is used to estimate the hydrophilicity of the 
facet. 


TABLE I: Sizes of the droplet and shearing simulations. Cell dimensions are measured in A. 


Interface 

Droplet 

^ice ^droplet 

Shearing 

^ice ^liquid 

Basal {0001} 

12960 2048 134.70 140.04 

900 1846 23.87 35.83 98.64 

Pyramidal {2021} 

11136 2048 143.75 121.41 

1216 2203 37.47 29.50 93.02 

Prismatic {1010} 

9900 2048 110.04 115.00 

3000 5464 35.95 35.65 205.77 

Secondary Prism {1120} 

11520 2048 146.72 124.48 

3840 8176 71.87 31.66 161.55 
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FIG. 2: Computational model of a slab of ice being sheared through liquid water. In this figure, the ice 
is presenting two copies of the prismatic {lOlO} facet towards the liquid phase. The RNEMD simulation 
exchanges both linear momentum (indicated with arrows) and kinetic energy between the central box and 
the box that spans the cell boundary. The system responds with weak thermal gradient and a velocity profile 
that shears the ice relative to the surrounding liquid. 
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TABLE II: Structural and dynamic properties of the interfaces of Ice-Ih with water. 


Interface 

Channel Size 

Width (A) Depth (A) 

Droplet 

9oo ( ) ^spread ) 

Shearing‘s 

^struct (-^) ^dyn (-^) 

Basal {0001} 

4.49 

1.30 

34.1(9) 

0.60(7) 

5.9(3) 6.5(8) 

3.2(4) 

2(1) 

Pyramidal {2021} 

8.65 

1.37 

35(3) 

0.7(1) 

5.8(4) 6.1(5) 

3.2(2) 

2.5(3) 

Prismatic {lOlO} 

6.35 

2.25 

45(3) 

0.75(9) 

3.0(2) 3.0(1) 

3.6(2) 

4(2) 

Secondary Prism {1120} 

6.35 

2.25 

43(2) 

0.69(3) 

3.5(1) 3.3(2) 

3.2(2) 

5(3) 


'^Liquid-solid friction coefficients (k^ and Kj,) are expressed in 10'^ amu A"^ fs h 
^Uncertainties in the last digits are given in parentheses. 
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I. THE ADVANCING CONTACT ANGLE 


The advancing contact angles for the liquid droplets were computed using inversion of Eq. (2) 
in the main text which requires finding the real roots of a fourth order polynomial, 

C4COS^0 + C3COS^0 + C2COS^0 + CiCOS0 + Co = 0 (1) 

where the coefficients of the polynomial are expressed in terms of the z coordinate of the center of 
mass of the liquid droplet relative to the solid surface, z = Zcm — ^surface, and a factor that depends 
on the initial droplet radius, k = The coefficients are simple functions of these two 

quantities. 


C4 = Z^ -h 

(2) 

C 3 = 8 z^ -f 

(3) 

C 2 = 24z^ + 

(4) 

Cl = 32 z3 

(5) 

cq = 16z^ -21 k^. 

(6) 


Solving for the values of the real roots of this polynomial (Eq. 1) give estimates of the advancing 
contact angle. The dynamics of this quantity for each of the four interfaces is shown in figure 1 
below. 

II, INTERFACIAL WIDTHS USING STRUCTURAL INFORMATION 

To determine the structural widths of the interfaces under shear, each of the systems was di¬ 
vided into 100 bins along the z-dimension, and the local tetrahedral order parameter (Eq. 5 in 
Reference 1) was time-averaged in each bin for the duration of the shearing simulation. The spa¬ 
tial dependence of this order parameter, q{z), is the tetrahedrality profile of the interface. The 
lower panels in figures 2-5 show tetrahedrality profiles (in circles) for each of the four interfaces. 
The q{z) function has a range of (0,1), where a value of unity indicates a perfectly tetrahedral 
environment. The q{z) for the bulk liquid was found to be ~ 0.77, while values of ~ 0.92 were 
more common in the ice. The tetrahedrality profiles were fit using a hyperbolic tangent function 
(see Eq. 6 in Reference 1) designed to smoothly fit the bulk to ice transition while accounting for 
the weak thermal gradient. In panels b and c of the same figures, the resulting thermal and ve¬ 
locity gradients from an imposed kinetic energy and momentum fluxes can be seen. The vertical 
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dotted lines traversing these figures indicate the midpoints of the interfaces as determined by the 
tetrahedrality profiles. 


III. INTERFACIAL WIDTHS USING DYNAMIC INFORMATION 

To determine the dynamic widths of the interfaces under shear, each of the systems was di¬ 
vided into bins along the z-dimension 3 A wide) and C 2 (z,t) was computed using only those 
molecules that were in the bin at the initial time. To compute these correlation functions, each 
of the 0.5 ns simulations was followed by a shorter 200 ps microcanonical (NVE) simulation in 
which the positions and orientations of every molecule in the system were recorded every 0.1 ps. 

The time-dependence was fit to a triexponential decay, with three time constants: x short, mea¬ 
suring the librational motion of the water molecules, Xmiddie, measuring the timescale for breaking 
and making of hydrogen bonds, and Xiong, corresponding to the translational motion of the water 
molecules. An additional constant was introduced in the fits to describe molecules in the crystal 
which do not experience long-time orientational decay. 

In Figures 6-9, the z-coordinate profiles for the three decay constants, Xshon, Xmiddie, and Xiong 
for the different interfaces are shown. (Figures 6 & 7 are new results, and Figures 8 & 9 are 
updated plots from Ref 1.) In the liquid regions of all four interfaces, we observe x^iddie and Xiong 
to have approximately consistent values of 3 — 6 ps and 30 — 40 ps, respectively. Both of these 
times increase in value approaching the interface. Approaching the interface, we also observe that 
X short decreases from its liquid-state value of 12 —16 fs. The approximate values for the decay 
constants and the trends approaching the interface match those reported previously for the basal 
and prismatic interfaces. 

We have estimated the dynamic interfacial width by fitting the profiles of all the three 
orientational time constants with an exponential decay to the bulk-liquid behavior, 

x{z) Xiiqnid{X-wall Xnqnid^S ^ wall)/ dyn 

where Xuquid and x^aii are the liquid and projected wall values of the decay constants, Zwaii is the 
location of the interface, as measured by the structural order parameter. These values are shown 
in table 1 in the main text. Because the bins must be quite wide to obtain reasonable profiles of 
C 2 (z,t), the error estimates for the dynamic widths of the interface are significantly larger than for 
the structural widths. However, all four interfaces exhibit dynamic widths that are significantly 
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below 1 nm, and are in reasonable agreement with the structural width above. 


[1] P. B. Louden and J. D. Gezelter, J. Chem. Phys. 139, 194710 (2013). 
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FIG. 1: The dynamic contact angle of a droplet after approaching each of the four ice facets. The decay to 
an equilibrium contact angle displays similar dynamics. Although all the surfaces are hydrophilic, the long¬ 
time behavior stabilizes to significantly flatter droplets for the basal and pyramidal facets. This suggests a 
difference in hydrophilicity for these facets compared with the two prismatic facets. 
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FIG. 2: Properties of the pyramidal interface being sheared through water at 3.8 ms"'. Lower panel: the 
local tetrahedral order parameter, q{z), (circles) and the hyperbolic tangent fit (turquoise line). Middle 
panel: the imposed thermal gradient required to maintain a fixed inferfacial femperafure of 225 K. Upper 
panel: fhe fransverse velocity gradienf fhaf develops in response fo an imposed momenfum flux. The verfical 
dolled lines indicate fhe localions of fhe midpoinls of fhe Iwo inlerfaces. 
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FIG. 3: The secondary prism interface with a shear rate of 3.5 ms"'. Panel descriptions match those in 
figure 2. 
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FIG. 4: The basal interface with a shear rate of 1.3 ms"'. Panel descriptions match those in figure 2. 
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FIG. 5: The prismatic interface with a shear rate of 2 ms'^ Panel descriptions match those in figure 2. 
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z-displacement from the center of the ice (A) 


FIG. 6: The three decay constants of the orientational time correlation function, C 2 {z,t), for water as 
a function of distance from the center of the ice slab. The vertical dashed line indicates the edge of the 
pyramidal ice slab determined by the local order tetrahedral parameter. The control (circles) and sheared 
(squares) simulations were fit using shifted-exponential decay (see Eq. 9 in Ref. 1). 
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FIG. 7: Decay constants for C 2 (z,t) at the secondary prism face. Panel descriptions match those in 6. 
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FIG. 8: Decay constants for C 2 {z,t) at the basal face. Panel descriptions match those in 6. 


12 

























FIG. 9: Decay constants for € 2 ( 1 , t) at the prismatic face. Panel descriptions match those in 6. 
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